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SOME PROBLEMS OF THE CALCULATION OF THREE-DIMENSIONAL 
BOUNDARY-LAYER FLOWS ON GENERAL CONFIGURATIONS 

by 

Tuncer Cebeci, Kalle Kaups, G. J. Mosinskis, J. A. Rehn 
Douglas Aircraft Company 

I . SUMMARY 

An accurate solution of the three-dimensional boundary-layer equations 
over general configurations such as those encountered in aircraft and space 
shuttle design requires a very efficient, fast, and accurate numerical method 
with suitable turbulence models for the Reynolds stresses. In the study 
reported here, we investigate the efficiency, speed, and accuracy of a numer- 
ical method together with the turbulence models for the Reynolds stresses. 

The. numerical method is an implicit two-point finite-difference method (Box 
Method) developed by H. B. Keller (ref. 1) and applied to the boundary-layer 
equations by Keller and Cebeci (refs. 2, 3). In addition, we study some of the 
problems that may arise in the solution of these equations. 

In Chapter II we write the governing boundary-layer equations, in both 
streamline and body coordinates, for three-dimensional, compressible, laminar 
and turbulent boundary layers. Those equations require initial conditions on 
two intersecting lines. Hence, we discuss the specifications of the initial 
conditions in longitudinal and lateral directions and the initial starting 
conditions, such as those existing in the nose region and those existing in 
the fuselage-wing junctures. We discuss the relative advantages of streamline 
and body coordinates and outline a solution procedure that combines both stream- 
line and body coordinates. 

When physical coordinates are used, the solutions of the governing boundary- 
layer equations are quite sensitive to the spacings in the streamwise direction 
(x), and the crosswise direction (z), and require a large number of x- and 
z-stations. In problems where computation time and storage become important, 
it is necessary to remove the sensitivity to ax- and Az-spacings. That can 
be done by expressing and by solving the governing equations in transformed 
coordinates. We, therefore, consider, in Chapter II, a convenient transforma- 
tion and express the boundary-layer equations in terms of transformed variables. 



In Chapter III we discuss Keller's Box Method. We investigate the computa- 
tion time and the accuracy of the method for two-dimensional and three- 
dimensional flows, and we compare the stability properties of the Box Method 
with the stability properties of the method used by Krause, Hirschel, and 
Bothmann (ref. 4). 

In Chapter IV we discuss and present a suitable turbulence model for 
three-dimensional compressible flows. The model, which is based on the con- 
cepts of eddy viscosity and eddy conductivity (turbulent Prandtl number), has 
given accurate results for two-dimensional incompressible and compressible 
flows and for three-dimensional incompressible flows. We also present results 
obtained with that formulation for an attachment-line, incompressible, turbu- 
lent flow on an infinite swept wing. 

Finally, in Chapter V we outline a procedure for computing the compres- 
sible three-dimensional multi-component-gas boundary layers on general 
configurations. On the basis of the studies conducted in the earlier chapters, 
we estimate the computation time and computer-storage requirements. 



II. GOVERNING EQUATIONS 


2.1 The Boundary-Layer Equations in Streamline Coordinates 

The governing boundary-layer equations for three-dimensional compressible 
flows in a curvilinear orthogonal coordinate system are given by the following 
equations: 

Continuity * 

§7 (ph 2 u) + (ph-jW) + |y (h^pv) = 0 , (2.1.1) 


x-Momentum 


u 3u , w 3u — — 3u 
B ST + p 57 + pv 57 


2 -Momentum 


puwK 2 + pw K-| 


L- i£ + 2- ( v 2“ 

h-j 3x 3y v 3y 


pu'v') ) 

( 2 . 1 . 2 ) 


U 3w , W 3W , — — 3w . 1 / , ,.2» _ 

0 h7 air + 0 pv 35-- puwK l + pU K 2 ■ 


L i£ + L ( u M 

h 2 3Z 3y ^ 3y 


pw'v') , 
(2.1.3) 


Energy 


u 3H , w 3H — 3H 
h^3T +pv 5y 


3_ 

sy 


u 

Pr 


3H , 1 \ 3 (u Z + vi 2 \ Trim 

37 + T ~ w) 37 \—r- ) ~ pv H 


(2.1.4) 


where pv = pv + p 'v* and h-j and h 2 are metric coefficients. The latter 
are functions of x and z, that is, 


h-| = h-j (x,z) , h 2 = h 2 (x,z) . (2.1.5a) 

The parameters K-j and K 2 are known as the geodesic curvatures of the curves 
x=const. and z=const., respectively. They are given by 


1 9h ? 

K-, = - ~L 

1 h-jl^ 3x 


1 3 ^1 
K 2 h-j h 2 3z 


(2.1.5b) 
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The streamline coordinate system is also an orthogonal coordinate system 
formed by the inviscid streamlines and their orthogonal trajectories on the 
surface. As is shown in figure 1, the projection of the free-stream velocity 
vector on the surface is aligned with the surface coordinate x. The velocity 
component along the z-axis, referred to as the cross-flow velocity is zero at 
the edge of the boundary layer. The x-momentum equation (2.1.2) is referred 
to as the streamwise momentum, and the z-momentum equation is referred to as 
the cross-flow momentum equation. 

At the edge of the boundary layer, (2.1.2) and (2.1.3) reduce to 


u 


s 




1 8p 
h-j 3X 


(2.1.6a) 


p s u s K 2 


1 3p 

h 2 


(2.1.6b) 
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The boundary conditions for the governing equations in streamline coordi- 
nates, (2.1.1) to (2.1.4), are 

y = 0 u,w = 0, v = v w (x,z) 

H - H w or (f£) - h; (given) (2.1.7a) 

y = 6 u = u s (x,z), w = 0, H = H g (2.1.7b) 

The solution of the system given by (2.1.1) to (2.1.4) requires closure 
assumptions for the Reynolds stresses appearing in these equations. They also 
require initial conditions on two intersecting lines. In some problems the 
initial conditions can be established with ease; in some problems they require 
careful studies. As an example, consider the blunted circular cone with 
an elliptic cap shown in figure 2. It is clear that before the points 
B 2 arrd C 2 can be calculated, it is necessary to calculate the initial 
points A-j , A 2 , and C-| . 



Figure 2. Blunted Circular Cone with Elliptic Nose Cap. Dashed Line 
Separates the Nose Cap from the Conical Configuration. 
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The initial points in the longitudinal direction, namely, A-j and 
can be calculated by taking advantage of the symmetry conditions since in that 
direction the flow is two dimensional except for the cross-flow derivatives. 

The flow in that direction is usually referred to as the attachment-line 
flow, because the attachment line is a streamline on the body on which both the 
cross-flow velocity components in the boundary layer and the cross-flow 
pressure gradient are identically zero. Since w, K 2 are zero on the attach- 
ment line, the cross-flow momentum equation is singular on that line. However, 
differentiation with respect to z will yield a nonsingular equation. After 
performing the necessary differentiation for the cross-flow momentum equation 
and taking advantage of symmetry conditions, we can write the governing 
attachment-line flow equations as: 

Continuity 


(ph 2 u) + ph.|W z + |- (h 1 h 2 pv) = 0 


Streamwise Momentum 


u 9u' — 3u u s 9u s A a / 9 u r\ 

p h, 9x p 9y p s h ] 9x 9y 9y pu v ) 


'1 


( 2 . 1 . 8 ) 


(2.1.9) 


Cross-Flow Momentum 
u 9w z A P 2 
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9K 


2 .2 


h, 9x + hT w z + pv 93r~ puK l w z = 9T- p s u s I 1 


pu 


2 

p u 
s s 


sy 


Energy 

u 9H , — — 9H _ 9 
p h 1 9x + pV 9y 9y 


9W 

t 

ay 


y VTT" ~ p(w* v 1 ). 


y 9H , /, 1 \ 9 /u 2 \ -rrrr 

Pr 9y + v \ ] Pr j 9y (2 j pv H 


( 2 . 1 . 10 ) 


( 2 . 1 . 11 ) 


where w z = 9w/9z. 

The attachment-line flow equations are subject to the following boundary 
conditions: 

y = 0 u,w z = 0 v = v w (x,z) 

H = H w or (§) = h ; (g^en) (2.1.12a) 
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y = s 


(2.1.12b) 


u = u g (x,z), w z = 0, H = H s 


The initial points A-| and A 2 can be obtained by solving the system 
(2.1.8) to (2.1.12). Again, however, the equations are singular at the 
first point (in our case A-| ) and require starting conditions. These and 
the initial conditions in the lateral direction, namely, B-j , , etc., can 

be obtained by constructing special solutions in the neighborhood of the 
stagnation point. At first, however, it is necessary to obtain the similarity 
equations known as the stagnation-point equations. They can be obtained by the 
procedure described below. 

The governing conservation equations for three-dimensional laminar 
compressible flows in rectangular coordinates are 
Continuity 


17 (pu) + gj (pw) + gy ( P v) = 0 


x-Momentum 


3 U , 3u , 3U 3p , 3 

pu 3 x t p, 82 + w 5y ° “)x )> 


(>») 


z-Momentum 


3w , . . 3w , 3w 3p , 3 / 3w \ 

pu 3i + pw 3? + pv ay = "3r + 37^ jy) 


Energy 


3H , ,, 3H , „ 3H 3 

pu 3x + pw 17 + pv 3y 5y 


2 . 2 


y 3H / n t \ 3 /u + w \ 
Pr 3y Pr/ 3 y \ 2 / 


We define a two-component vector potential by 


pu “ I? ’ 


pW 


_ 3 if) 

~ 3y ' 


P v 


= _ (it + i±\ 

\3X 3Z/ 


and introduce the following transformation 


(2.1.13) 


(2.1.14) 


(2.1.15) 


(2J.16) 


(2.1.17a) 


d « = p o p o u e dx ’ 


dn - 


pu. 


(2c) 


Tff 


dy 


(2.1.17b) 


together with two scalar functions ip and <p defined by 
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(2.1.17c) 


♦ = (2 5 ) 1/2 f( n ), 4 . = <2S) 1/2 g(n) 

e 

The subscript 0 denotes the reference conditions, and the subscript e 
denotes the edge conditions. 

Introducing the relations given by (2.1.17) into (2.1.14) to (2.1.16), 
we get 


(Cf")' + p 6 /p - (f 1 ) 2 + ff" + cgf 11 = 0 (2.1.18) 

(Cg")' + c[p e /p - (g') 2 + gg"] + fg"= 0 (2.1.19) 


(fe 0 ')' + (f + Cg)6, = ° (2.1.20) 

where the primes denote differentiation with respect to the similarity 
parameter n. Those equations assume that the outer flow is irrotational and 
that its components, upon suitable rotation of coordinates, are given by 
u g = Ax and w g = Bz, where the constants A and B are related to the 
shape of the body near the stagnation point. In addition, the parameters 
C, c, f', g' and 0 are defined by the following expressions: 


c = ft- • c = B/A ’ f ' - s- • 9' - ir • 0 - H/H e 

M o M o e e 


( 2 . 1 . 21 ) 


The system given by (2.1.18) to (2.1.20) is subject to the following 
boundary conditions: 


n = 0 f = f* = g = g' = 0, 0 = 0 W or 0^ = given (2.1.22a) 

n = noo f* = g’ = 1 0=1 (2.1.22b) 

Once the system (2.1.18) to (2.1.22) is solved, the initial conditions in 

the streamwise and the crosswise directions can be obtained in the following way 

Let us use x and z to denote the location at which we want to specify 
00 

the profiles. The velocity profile making an angle a with the x-axis is 


u = u f' cos a + w e g' sin a 
and the external flow component u g is 
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u s = 

u e 

cos a + w„ sin a 
e 


Then the nondimensional 

streamwise 

profile i*s 



u 

f' 

+ c(z 0 /x 0 )g' tan 

a 


u s " 

1 

+ c z 0 /x 0 tan a 


If a is associated wi 

th the 

streamline direction 

(see figure 

expression becomes 






u 

f 

’ + c2(z o /x o )2g ’ 



u s 


1 ♦ c 2 (z 0 /x 0 ) 2 


Similarly, the cross-flow 

f component 

w at the same 

location is 


w 

c 

V X o (g ’ - f '> 



u 

s 

1 

+ c 2 (z 0 /x 0 ) 2 



(2.1.23) 


(2.1.24) 



Figure 3. Resolution of the Velocity Profiles Near the Stagnation Point into 

Streamwise and Cross-Wise Components, tan a = w /u =c z /x . 

, e e oo 
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The velocity profiles given by (2.1.23) and (2.1.24) can now be used as 
initial velocity profiles at . They can also be used to represent the • 
initial profiles at (note that w = 0 now). However, better initial 
profiles in the neighborhood of the stagnation point can also be obtained by 
following the procedure discussed by Squire (ref. 5). Once the profiles at 
A-j are known, the attachment-line flow equations (2.1.8) to (2.1.12) can be 
solved to obtain the solution at A^. Then it is obvious that the general 
streamline equations (2.1.1) to (2.1.4) can be solved subject to the boundary 
conditions (2.1.7). 

There are other practical problems in which the initial conditions as 
described in figure 1 cannot be obtained as readily. As an example, consider 
figure 4. Here, the initial conditions require special considerations. 
Clearly, the attachment line AB along the wing leading edge is a plane of 
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symmetry, and the initial conditions on AB can be calculated as before by 
solving the attachment-line flow equations. However, the initial conditions 
on AC that form the wing-fuselage juncture cannot be calculated easily. 

In fact, the viscous flow along the line AC is not of the boundary-layer 
type. It belongs to a class known as the boundary-region type. Certain 
approximations must be made to specify initial conditions on that line. 

2.2 The Boundary-Layer Equations in Body Coordinates 

2,2.1 Remarks on the Two Coordinate Systems 

The discussion in Section 2.1 concerning figure 1 was based on the 
assumption that calculations from the initial data lines are to proceed in a 
streamline coordinate system. Although the streamline coordinate system is 
very general , its calculation is a major undertaking in itself and must be 
repeated at every change of attitude of the body. If the body is relatively 
simple, advantage can be taken of a geometry-oriented (body) coordinate 
system, which will eliminate the need to calculate the streamlines for 
each change of attitude or Mach number. The only disadvantage of body 
coordinates is that the initial lines cannot always be made to coincide with 
the body coordinate lines. For example, on a body of revolution (see fig. 5) 



Figure 5. A Body of Revolution at an Angle of Attack 


n 


at an angle of attack, the stagnation point S is removed from the nose, 
which is the origin of the body-coordinate system. In order to start the solu- 
tion, it is necessary to calculate the stagnation-point flow along the line 
BC. To proceed further, it is convenient to use streamline coordinates for a 
short distance in order to avoid the body singularity at A. The change from 
streamline coordinates to body coordinates should be done at a constant value 
of one body coordinate, for example, on line ED. The location of the line 
ED is arbitrary as long as point D is aft of point C. As shown in figure 5, 
the inviscid streamlines do not necessarily intersect the line ED at a con- 
stant body-coordinate interval. For that reason, interpolation of the rotated 
profiles (to be discussed in Section 2.2.2) on the line ED is unavoidable if 
calculations with constant increments in the body coordinates are desired. 


2.2.2 Equations in Body Coordinates 

The governing boundary-layer equations for three-dimensional compressible 
flows in body coordinates are given by (2.1.1) to (2.1.5). 

At the edge of the boundary layer, (2.1 .2) and (2.1.3) reduce to 


u„ 

e e 

5 e h-j 3x 


w e 8u e 2 1 9p 

’e 3T~ - °e u e w e K 2 + p e w e K l = " 


(2.2.1a) 


U„ 3W W 3W 

P tr — tt — — + p 7-— - — — - p U W K, + p U = 
e h-j 3x H e h^ az e e e 1 e e 2 


1 3p 

h 2 3z 


(2.2.1b) 


The boundary conditions for the governing equations in body coordinates are 
y = 0 u,w = 0, v = v w (x,z) 

H = V or (!y) w = H w (given) (2.2.2a) 

y = <5 U = u w (x,z), w = w e (x,z), H = H e (2.2.2b) 

Making use of the symmetry conditions, we can write the two attachment- 
line flow momentum equations as 
x-Momentum 



z-Momentum 


u 9w z 12 — 9w z 2 9K 2 

p h7 *T + p h7 w z + pv 8F “ puK l w z + pu JT 


'1 


1 3 P 

ho 


'2 9z 


3y 


9W 

" ~ P(v7 T V r ) z 


(2.2.4) 


The continuity and the energy equations are still the same, (2.1.8) and 

(2.1.11) , respectively. Similarly, the boundary conditions are the same as 

(2.1.12) , except that now the subscript s on u and H should be replaced 
by e. 

The solution of the governing equations in body-coordinates aft of line ED 
(see fig. 5) requires initial velocity profiles, which come from the solution of 
the governing equations in the streamline coordinates. Except for the attach- 
ment line, they can be obtained in the following way. 

Let us write the velocity components in streamline coordinates with 
bars, namely, u", w, and the angle the external streamline makes with the body 
coordinate x-direction as y (see fig. 6). Then the velocity components u 
and w in the body-coordinate system x and z are 

u = lT cos y — w sin y (2.2.5a) 

w = IT sin y + w cos y (2.2,5b) 



Figure 6. Notation for Streamline (Barred) and Body (Unbarred) Coordinates. 
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where 


. -l/ w e 
Y - tan f- 

\ e, 


( 2 . 2 . 6 ) 


On the attachment line the coordinate directions coincide, so that 

u = U and w = w = 0 


An expression for w z can be obtained by making use of the expressions 


w 

sin Y = / , 
s 


u 

cos Y = ^ 


and by taking limit of (2.2.5b) as z -> 0. The result is 


_ aw_ _ u_ 9w e 8w dz 

W z = 3Z " u s 32 sT dZ 


But, in the limit. 


or 


h^dz .= I^dz 


W 

z u c 3Z - t- 
s 9 z 


(2.2.7) 


2.3 Transformation of the Governing Equations 

In this, section we shall consider the transformed form of the governing 
equations discussed in the previous two sections. Although those equations 
can be solved in their physical coordinates x, y, z, it is often convenient 
to solve them after they have been expressed in terms of transformed coordinates. 
In problems where the computer storage becomes important, the choice of using 
transformed coordinates becomes necessary, as well as convenient, since the 
transformed coordinates allow large steps to be taken in the x and z direc- 
tions. The reason is that the profiles expressed in the transformed coordinates 
do not change as rapidly as they do when they are expressed in physical coordi- 
nates. In addition, the use of transformed coordinates stretches the coordinate 
normal to the flow and takes out much of the variation in boundary-layer thick- 
ness for laminar flows. 
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2.3.1 Streamline Coordinates 

Following Moore (ref. 6), we define a two-component vector potential such 


ph 2 U = fy* ph l W “ ly (2.3.1a) 

h 1 h 2 PV = “(ff + |i) + h l h 2 (pv) w (2.3.1b) 


We note that equation (2.3.1b) includes the effect of mass transfer and 
decouples the wall boundary conditions. 

We also define the following transformations: 

/ u s \ 1/2 

x = x, Z = z, dn = ( p -~— ) pdy (2.3.2a) 

4> = ( p s M s u s x ^^ h 2 f(x, z, n) (2.3.2b) 

<P = (p s p s u s x ) 1/2 h ig (x, z, n ) (2.3.2c) 

Introducing the expressions given by (2.3.1) and (2.3.2) into (2.1.2) to (2.1.4) 
and making use of the relations given by (2.1.6), we get 

Streamwise Momentum 



(2.3.3) 

Cross- Flow Momentum 




Energy 
C 


u. 


1 


+ P 0 fe' + (R + N + 2M) 


(PV) W r 1/2 , _ X ( f , 86 , 3fV X / , 86 , S£^ 

s ~S\ ** **/ h 2\ ^ 


ge ' 
2hT 


(2.3.5) 


where the primes denote differentiation with respect to n and 


f 1 = — 
r u » 
s 

3 U 

P = — S 


1 - U s 3x • 


N = 


x_^s_ 
u$ 3 z 


i w 

« = - r • 

s 


0 - iL 

9 - H s * 


u s x 

R = — , 
s v s 


C = PM 


p s y s 


M H - K 2 h 2* S = R = (p s U s ) 


S"S 


P 2 ^ + P 1 + S 2K i h i x ) 2h, 


(2.3.6) 


In the above equations we have used Boussinesq's eddy-viscosity and eddy- 
conductivity concepts in order to satisfy the closure conditions for the 
Reynolds stresses. They are defined by 


-pu v 


pe 


3U_ 
x 3 y 


-pw 1 v 1 = 


pe T 


3 W 

ay 


- p H ' V 1 = pe 


3H 

0 ay 


(2.3.7a) 


The turbulent Prandtl number and the dimensionless transport coefficients are 
defined by 


+ 




e 0 



V 


(2.3.7b) 


The boundary conditions (2.1.7) become 


n = 0 

f = g = 0, 

f = g* = o. 

CD 

II 

CD 

Z 

or e' = e' 
w 


f' + 1 

g' o 

e - l 

(2.3.8) 


The attachment-line equations can also be transformed by a similar 
procedure. This time, we define the two-component vector potential by 


ph 2 u 



ph l w z 


9<j> 

ay ’ 


h-jl^pv = — 




+ h 1 h 2 ( P v) w 


(2.3.9) 


16 



and again use the expressions given by (2.3.2). Introducing the expressions 
(2.3.9) and (2.3.2) into (2.1.9) to (2.1.11), we get 
Streamwise Momentum 


[C(l + + J- 


P, ? 

r- (f,) 


+ P,ff" + gf" - 


!pv) 


w pi /2^„ 


p s u s s 


1 


Cross-Flow Momentum 


V 3f 

, ax ax 


(2„3. 10) 


[C(l + e +)g"]' + ^ [gg“ - (g 1 ) 2 ] + P 2 fg” + (y 

+ [^-(f) Z l-^Ry 2 g" -i-(* 


» rV 2 „ . x L, aal „ il 


p s u s s 


ax J ax 

/ 

(2.3.11) 


Energy 
C 


1 E Pr. JPr 6 H 

< t / s 




(pv) 


+ p 2 fe ' + ^ge' 


p s u s s 


w r 1/2 0 . = X f , M_ 0 , |f 


ax 


ax 


(2.3.12) 


where the definitions of the terms are the same as those defined in (2.3.6), 
except for g', which is equal to w z /u s . 

The boundary conditions (2.1.12) become 

n = 0, f = g = 0, f* = g' = 0, 0 = e w or 0 1 = 0 ^ (2.3.13a) 

n - f' = e = 1 , g ' = 0 (2.3.13b) 

2.3.2 Body Coordinates 

The relations used to transform the equations in body coordinates are 
similar to those used in the previous section. For the general case, we again 
use the two-component vector potential defined by (2.3.1) and the same rela- 
tions defined by (2.3.2a,b), except that now the subscript s is replaced by e, 
that is. 
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(2.3.14) 


dn = fevT P dy 


1 /2 

* = (p e v e u e x) h 2 f(x, z, n) 


and (() is defined by 


, 1/2 


r w 


= (p e y e u e x ^ h l liT ) 9 ^ Xj z ’ 


( 2 .: 


( 2 . 


With these relations and with those given by (2.2.!), we can write (; 
to (2.1.4) as 
x-Momentum 


[C(l + ejjf"]' + P 2 ff" + P 3 gf" + l 
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z-Momentum 
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Energy 
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3.15b) 
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where 


f* = u/u , g' = w/w . e = H/H , R = u x/v 


e 

9U „ 3l L 

M = — , N = e 


3x 

e 


U 9 Z ’ 

e 


e’ x 
aw 

p = x e 

w 9Z 

e 


e e 


Q = 


3W 

x e 

w„ 9x 
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S ~ p e y e 3x 




P 2 = (1 + M + S - 2K- | h. | x) 


( 2 . 


P = 1 ( 2p „ N + R _ 2K 9 h,x) 

3 u e 2h 2 2 L 


v u e/ 1 


P 5 = ” u~ h ? ( “ K 2 h 2 x + N) ’ 
e 2 


P c = — K 9 x 

6 w e 2 


P 7 (K^h^x Q) I 


1 

The boundary conditions (2.2.2) become 

n=0 f = g = 0 f 1 = g' = 0 

n n f 1 = g 1 = 1 0 = 1 


0=0 or 0 ' = e' (2. 

w w 

( 2 . 


The attachment-line equations can also be transformed by a similar 
procedure. We define the two-component vector by the relations given by 
(2.3.9) and again use the relations (2.3.14), except that now we 
def i ne <j> by 

1 /2 W 7 

4 > = (P e p e u e x ) 7 h l 7^ g(x ’ z ’ ^ 

e 

Introducing the expressions (2.3.9), (2. 3. 14a, b), and (2.3.21) into (2. 
(2.1.27), and (2.1.11), we get 
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x-Momentum 
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z-Momentum 


[c(l + 4>g"]' + P 2 fg“ + jp gg" + p 8 (f'g' 

+ p g[r- (f,)2 ]-57 R x /2f " 


^^f^-(g') 2 


h 2 p 


p u x 
e e 


w R l/2 flI x fl 3^_ g „ 


Energy 


C 1 + 


Pr t/ Pr H e ' P r) 


x / f) |l_ e . 


(2.3.23) 


(pv) w 1/2 

+ P 2 fe ' +H -ge' -p^R^e' 

^ e e 

(2.3.24) 


f' = u/u , 


g' = w z /w ze , 


e = H/H 


P 8 = ( K l h l x w ze af)^ 


3 W 

p _ x e 

*1 “ u~ W 
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xu 3 1C 

p 9 = i^?r < 2 - 3 - 25 > 

ze 


The boundary conditions are 
n = 0 f = g = 0 

n = n f' = g' = 0 = 1 


f ' = g ' = 0 0 = 0 or 0' = e, , 

3 w w 


(2.3.26b) 
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III. KELLER'S BOX METHOD 


The governing boundary-layer equations presented in the previous chapter 
form a system of coupled nonlinear partial differential equations that are quite 
difficult to solve. For a multicomponent gas, their solution is even more dif 
difficult because, in addition to mass-continuity, momentum, and energy equa- 
tions, we have a number of species-continuity equations to consider. The 
solution of those equations for general configurations such as those that^ occur 
in aircraft and space shuttle design requires a very efficient, fast, and 
accurate numerical method with suitable models for the Reynolds stresses. 

In this chapter we shall discuss the efficiency, speed and accuracy of a 
two-point finite-difference method developed by H. B. Keller (ref. 1) and applied 
to the boundary-layer equations by Keller and Cebeci (refs. 2,3). We shall 
investigate the computation time and accuracy and the stability properties of 
this method for two-dimensional incompressible laminar and turbulent flows 
as well as for three-dimensional laminar flows. On the basis of that informa- 
tion, in Chapter V, we shall estimate the computation time for three-dimensional 
laminar and turbulent boundary layers of a multicomponent gas, and we shall out- 
line an efficient procedure for solving those equations. But, first, we shall 
present a brief description of the Box Method and point out the several advan- 
tages of that method over the numerical methods now being used for boundary-layer 
calculations. For simplicity, we shall consider the infinite swept-wing equa- 
tions for an incompressible flow. 

3.1 Box Method for Infinite Swept-Wing Equations 

The transformed boundary-layer equations for an incompressible flow over 
an infinite swept wing follow from (2.3.16) and (2.3.17). With h-j = h 2 = 1 
and spanwise derivatives of the form d/dz being zero, they are 

Chordwise Momentum 0 I i-f • \ 

(bf")' + P 2 ff" + M[1 - (f'r] = x (f 1 j~~f n ~) (3.1.D 

Spanwise Momentum 

(eg-)' + ^ fg" = *(r If -g“ f) <3.i.2) 

where 

b = 1 + £* , c = l + e z , P 2 = ~ ~~2~ ~ (3.1.3) 
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We first write the two momentum equations in terms of a first-order 
system of partial differential equations. For that purpose we introduce new 
independent variables u(x,n)j v(x,n), w(x,n)» and t(x,n) so that we can 
write (3.1.1) and (3.1.2) as 


f 1 = u (3.1 .4a) 

u' = v (3.1 ,4b) 

g 1 = w (3.1.4c) 

w' = t (3.1. 4d ) 

(bv)' + P 2 fv + M(1 - u 2 ) = x (u |j- v |“ j (3.1 .4e) 

(ct) ■ +7 ft = x(u f£-tf±) (3.1.4f) 


We next consider the net rectangle shown in figure 7. We denote the net 
points by 


X 

o 

II 

o 

V# 

x n = Vl + k n’ 

n = 1, 2, 

..., N 

(3.1.5) 

n Q = 0 

n j = n j-l + h j’ 

j = 1, 2, 

. . . i J 

= n co 



Figure 7. Net Rectangle for the Difference Equations 
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The net spacings, k and h-, are completely arbitrary and indeed may have 

* ' vJ 

large variations in practical calculations. Such flexibility is especially 
convenient in turbulent boundary-layer calculations, which are characterized by 
large boundary-layer thicknesses. To get accuracy near the wall, small net 
spacing is required; large spacing can be used away from the wall. 

We approximate the quantities (f, u, v, g, w, t) at points (x n , n ^ ) of 
the net by net functions denoted by (f^, u^, v'?, g 1 ?, w?, t^). We also employ 
the notation, for points and quantities midway between net points and for any 
net function q*J: 

x n-l/2 - I ( s n + S n~l ^ ’ n j-l/2 " 7 + n j-l ) 

(3.1.6) 


n-1/2 _ 1 f n . n-1 n 

q j = ? (q j + q j >• 


n 1 / n , n \ 

q j-)/2=2 (q j +q j-l> 


The difference equations that are to approximate (3.1.4) are now easily formu- 
lated by considering one mesh rectangle as in figure 7. We approximate 
(3.1.4a) to (3.1 ,4d) using centered difference quotients and average about the 


midpoint (x 


n* 


n 0-l/2 ) 


of the segment P 2 P 4 , with the following results: 


f n 


f" i 
J-l _ 



h. 


J 
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u . 

— u . 
J-l 
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h. 


J 

n 

- n 

9 j 

“ 9 J-l 


h. 


J 

n 

n 

w . 


J 


= u 


h o 


j-l/2 


n 

V j -1/2 


W j - 1 / 2 


t n 

J-l/2 


(3.1.7a) 


(3.1.7b) 


(3.1.7c) 


(3.1 .7d) 


Similarly (3.1.4e,f) are approximated by centering on the midpoint 
x n-l/2’ n j -1/2 of ttie rectangle p -j P 2 P 3 P 4 ’ which gives 
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(bv)7 - (bv)" , - „ 

— Sr + •n»“ 2 )J. 1 / Z + + «„Kfv)".i / z 

, n-1 , n .n-1 x _ T n-1 

a n (f j-1/2 v j-l/2 v j-l/2 f j-1/2 - T j —1/2 


(3.1 ,7e) 


(Ct)" - (ct)" , , 

J 1 + (I + a)(ft)" 


S' 


j-1/2 “n (uw) j-l/2 


/ ,.n-l n . ,n-l n .n-1 .n 

a n U j-l/2 W j-l/2 W j-l/2 U j-l/2 t j-l/2 f j-l/2 


+ f n “^ t n ) = S n_1 f3 1 7fl 

j-1/2 j-l/2 ; j-1/2 


where 


T n_1 

'j-1/2 


(f v ) n "l _ ( u 2 l n_1 
™ j-1/2 lU j-1/2 


M n 

— M — 
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h j 


, M n-1 I , /, ,2,n-l . D n-1 ,n-l 

+ H l’ ~ (u ) j-1/2 + P 2 fv j-1/2 


^n-1 

j-1/2 


= a 


( ft )j-l/2 (uw) j-1/2 


n- 1/2 


(3.1.8a) 


x —— x , 
n n- 1 


(ct )"- 1 - (ct)"-; , 

^ FT ^ + 2 ( ft ,j-l /2 j 


(3.1.8b) 

(3.1.8c) 


Equations (3.1.7) are imposed for j =1,2, ...» J. For most laminar flows 
Dj is constant. For turbulent flows, nj may be increased, with no essential 
difficulty, as the calculations proceed downstream from the point of transition. 
The boundary conditions for (3.1.1) and (3.1.2) are 






u 0 ' 0- 


w „" ■ 0. 


I ' 1 

Uj = 1, 


9? = 1 


(3.1.9) 
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If we assume (f?~\ u!?"\ v 1 ? - "^, g?~^ , w i~^ » ^i"^) t0 known for 
0 <_ j ^ J, then (3.1.7) for 1 < j < J, and the boundary conditions (3.1.9) 
yield an implicit nonlinear algebraic system of 6J + 6 equations in as many 
unknowns ( f !? » u" v?, g?, w?, t?). The system can be solved very effectively 

J J J J J J 

by using Newton's method. The details are presented in reference 3. The 
important observations are that the linearized equations obtained by applying 
Newton's method to (3.1.7) and (3.1.9) form a block tridiagonal system (with 
6x6 blocks) and that system can be solved very efficiently by the procedure 
discussed in reference 3. 

3.2 Computation Time of the Box Method 

We have studied the computation time of the Box Method for two-dimensional 
laminar and turbulent flows as well as for three-dimensional laminar flows. 
These studies were made on an IBM 370/165. 

From a computational aspect, turbulent boundary layers present a much more 
difficult problem of calculation than laminar boundary layers. Consider, for 
example, an incompressible turbulent flow. The skin-friction is appreciably 
greater than it is for a laminar flow yet the boundary-layer is much thicker. 
This means that the velocity gradient au/gy is greater at the wall. To 
maintain computational accuracy when au/ay is large, short steps in y must 
be taken; when it is small, longer steps can be taken. Therefore, near the 
wall the steps in a turbulent boundary layer must be shorter than they. are in 
a laminar boundary layer under similar conditions, yet near the outer edge 
they can be longer. 

The numerical method described in Section 3.1 is unique in that various 
types of spacings in both x- and y-directions can be used with ease. In 
the calculations we present in this chapter, we have done the calculations for 
an arbitrary Ax-spacing but for a particular An-spacing. The net in the 
n-direction is a geometric progression having the property that the ratio of 
lengths of any two adjacent intervals is a constant; that is, h. = Kh. , . 

J J 

The distance to the j-th n-line is given by the following formula: 

n j = h l K-J J = 1, 2, 3, ... J, K 2 < 1 (3.2.1) 
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There are two parameters: h-j , the length of the first Ao-step, and K, the 

ratio of two successive steps. The total number of points J can be calculated 
by the following formula: 


J = 


An[l + (K - 1) n 
£n K 


(3.2.2) 


In our calculations we select the parameters h-| and K and calculate 
the n . 

oo 

To study the computation time of the Box Method for two-dimensional 
turbulent boundary layers, we selected a flat-plate flow. In the range of 
Reynolds number R between 1 x 10 to 40 x 10 , 21 x-stations and 50 
n-points were computed. The total Central-Processing-Unit time (CPU) was 
0.048 min. That time corresponds approximately to 0.14 sec/x-station and to 
2.75 x 10“ 3 sec per n-point and per x-station. In the calculations the wall 
shear parameter, f^, was taken as the convergence parameter. The itera- 
tions were repeated until 


f ..(v+l ) _ f „(v) 

w w 

i /VrV M ( ^ + -p" 

' L w w 


] 


Y i 


(3.2.3) 


where y-j is a small error tolerance parameter. On the average, the calculations 
required two iterations per x-station with y^ - 0.01. 

To study the computation time of the box scheme for a semi-three-dimensional 
flow, we have considered two different test cases. In one test case we have 
computed the turbulent boundary layers over a yawed flat-plate approximately 
in the same Reynolds number range as the two-dimensional test case discussed 
above. Again 21 x-stations and 50 n-points were taken. The CPU time was 

0.085 min. That time corresponds approximately to 0.243 sec/x-station and 

-3 

4.8 x 10 sec per n-point and per x-station. Again in each x-station, calcu- 
lations required two iterations to satisfy (3.2.3). 

In the second test case we have considered the Bradshaw-Terrell flow 
(ref. 7), which is a flow past a 45° "infinite" swept wing. In- that flow, 
measurements were made only at the rear of the wing in a region of nominally 
zero-pressure gradient and decaying crossflow. The CPU computation time to 
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calculate the complete flow field with 20 x-stations and 50 n-points was 

0.051 min. That time corresponds to approximately 0.14 sec/x-station and 

-3 

3 x 10 sec per n-point and per x-station. At first, the computation time 
for this flow appears to be approximately the same as the time for the two- 
dimensional flow which required 0.048 min. However, the Bradshaw-Terrell 
flow required approximately one iteration per x-station and 22 iterations for 
the complete flow. Figure 8 shows a comparison of calculated and experimental 
results. 

The computation time of the Box Method was also studied for a full 
incompressible three-dimensional laminar flow by considering the flow past a 
flat plate with attached cylinder (see fig. 9). In this case, the iterations 
were repeated until 


fll(v + l ) JMl(v) 

r w r w 


0.0001 


(3.2.4) 



Figure 8. Results for the Relaxing Flow of Bradshaw and Terrell . The Calcula- 
tions Used the Eddy Viscosity Formulation Described in Chapter IV. 
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z 


Figure 9. Flow Past a Flat Plate with Attached Cylinder. 

The total CPU time for 25 x-stations, 16 z-stations and 21 n-points was 1.285 

_3 

min. That time corresponds to approximately 3 sec/plane and 9 x 10 sec per 
x-station, per z-station, and per n-point. On the average the calculations 
required 2 to 3 iterations on the attachment line and only 2 iterations away 
from the attachment line. The results agree quite well with those obtained 
by Dwyer (ref. 8) and by Fillo and Burbank (ref. 9), who have also studied 
the same flow using a different finite-difference method. 

3.3 Accuracy of the Box Method 

The accuracy of the Box Method has been studied for both incompressible 
and compressible, laminar and turbulent boundary layers past two-dimensional 
and axi symmetric bodies. Some of the results have been reported in refer- 
ences 2 and 3 and others will be reported in a forthcoming book by Cebeci 
and Smith. The results indicate that the method is quite accurate and 
extremely well suited for boundary layers, especially for turbulent flows. 
Extensive studies with incompressible and compressible turbulent boundary 
layers show that, in general, 40 to 50 n-points with the Box Method give 
results which are comparable to the results obtained by the method of 
reference 11, using 300 to 400 n-points. 

The studies in two-dimensional flows also show that one can take relatively 
large 4x-spacing in the x-direction as long as the equations are solved in terms 
of the similarity variables similar to the ones discussed in Chapter II. In 
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general, an airfoil calculation in transformed coordinates requires 20 to 25 
x-stations. However, the same calculation in physical coordinates may require 
50 to 75 x-stations. 

To study the effect of An- and Ax-spacings on the results, we have 

computed the turbulent boundary-layer flow on a flat plate for a Reynolds 

6 9 

number range of 10 to 10 . Figure 10 and Table 1 show the skin-friction 
results with two different An- and Ax-spacings. Figure 10 shows the results 
with fixed An-spacing (h, = 0.002, K = 1.226), and with variable Ax-spacing. 

* fT 

The latter was chosen such that starting from R = 10 , the aR -spacing of 

A X 


4-i 


° AR X = 2 2 X 10 
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Figure 10. Effect of Ax-Spacing on the Computed c -Results. Calculations Were 
Made for a Fixed An-Spacing. ' 


Table 1. Effect of An-Spacing on the Computed Results with a Fixed Ax-Spacing. 
h{ 0) = 0.002, h^ 0) = 0.001. 


R )0' 6 

X 

1 

c f (h{ 0) )10 3 

c f (h 1 ( 1 } }10 3 

R (hj°^)10 3 

Mhj^Jio 3 


R (h|°)h| 3 ^ ) 1 0 3 

1 .0 

3.583 

3.570 

2.23 

2.22 

3.566 

2.2167 

10.7 

2.387 

2.369 

15.2 

15.1 

2.363 

15.0667 

115.3 

1 .745 

1.731 

115.9 

115.1 

1 .726 

114.83 

1133.5 

1 .352 

1.329 

864.0 

850.9 

1 .321 

846.53 
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2 n ^ 2 x TO 6 and 2 n ^ x 10 6 gives approximately 20 and 40 x-stations, respect- 
ively, in the Reynolds number range under consideration. The results indicate 
that the c ^-values are not very sensitive to the Ax-spacing. 

Table 1 shows the computed and R„ values for fixed Ax-spacing 
n/4 6 t 0 

(aR = 2 ' x 10 ) with variable An-spacing. The calculations were first 

X 

made with h-j = 0.002, K = 1.226 and then the net points in the n-direction 

were doubled by halving each An-interval. 

Table 1 also shows the Richardson-extrapolated values of c* and R_. 

T (0) 0 

According to the results, the c^ and R 0 values computed by hj = 0.002 
spacing (approximately 50 n-points across the boundary layer) are quite 
satisfactory. 

Figure 11 shows the computed transformed bound ary- layer thickness, n^. as 
the calculations proceed downstream (see ref. 10). It is interesting to note 
that although the increases from 20 to 200, the use of the variable grid 
keeps the number of n-points approximately constant, and the use of the Box 
Method maintains the computation accuracy in a large range of Reynolds numbers. 

3.4 Stability Properties of the Box Method 

Currently there are a number of numerical methods used to solve the 
boundary-layer equations. In reference 11, Blottner gives a general review 
of these methods. The stability properties of most of these methods, except 



Figure 11. Variation of the Transformed Boundary-Layer Thickness and the 
Number of n-Points with Reynolds Number. 



for the method of Krause, Hirschel and Bothmann (ref. 4), have not been 
investigated for three-dimensional boundary-layer equations. 

It is difficult to compare the stability properties of the schemes of 
Krause, et al , with those of the box-scheme by the methods used in reference 4. 
One main reason is that Krause, et al, do not give a complete stability analysis 
but study only one momentum equation while neglecting some of the convection 
terms and other coupling terms. The box scheme is based on a different formu- 
lation of the boundary-layer equations (requiring them to be replaced by a 
first-order system and using transformed variables to reduce the variations in 
the solutions). However, when an analogous linearized stability analysis is 
made of the box-scheme (ref. 12), dropping terms similar to those neglected in 
reference 4, the stability properties are found to be at least as good as those 
of the best scheme of Krause, et al. Indeed even retaining terms that are 
dropped in reference 4, the analysis shows stability under very general 
conditions. 

Unfortunately, such linearized stability studies cannot be conclusive. 

The best test would be to solve several difficult problems with each method. 

We could not make such a comparison in our study because neither the exact 
difference scheme nor the exact problems treated were specified in reference 4. 

However, some important comparisons can be made. Since Krause, et al, 
always employ three-point differences in the normal (or boundary layer) direc- 
tion, they must use a uniform net through the boundary layer or else they do 
not get second-order accuracy. The box-scheme is unrestricted in net spacing, 
getting not only second-order accuracy but even fourth-order or sixth-order 
accuracy with only one or two Richardson extrapolations, respectively (ref. 2,3). 
Also, the Newton iterates used to solve the nonlinear (implicit) equations of 
the box-scheme converge quadratical ly and thus are very efficient and do not 
degrade the accuracy of the solution. It is never clearly spelled out how the 
nonlinearities are treated in reference 4, so comparisons here are again 
difficult. Finally, we note that the most stable scheme described by Krause 
et al does not have second-order accuracy in both tangential directions unless 
the net spacing is uniform in an appropriate one of these tangential coordinates. 
Again, there is not such restriction on the box-scheme. 

A complete analysis of the three-dimensional boundary-layer equations has 
never been made. But a preliminary investigation indicates that they are not 
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stable or well-posed for all tangential flow fields. Indeed something like 
this must be true since even in two-dimensional flows the boundary-layer 
equations become unstable when the tangential velocity changes sign (i.e., at 
separation). For flows in which the boundary- layer equations are not 
well-posed, it is impossible to devise stable and accurate difference schemes. 
If the tangential component of the velocity vector turns through a sufficiently 
large angle, this phenomenon seems to occur. This question should be studied 
in more detail in order to devise numerical schemes of maximum stability, or 
indeed to verify if the box-scheme, or any other scheme, possesses maximum 
stability properties (i.e., is stable whenever the boundary-layer problem is 
well posed). 
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IV. TURBULENCE SHEAR MODELS FOR THREE-DIMENSIONAL BOUNDARY LAYERS 


The solution of the boundary-layer equations for a turbulent flow requires 
closure assumptions for the Reynolds stresses. That can be done by a number of 
approaches*. One approach is to use simple eddy-viscosity and mixing-length 
formulas for the Reynolds stresses. The methods that use that approach are 
called mean-flow methods. Typical examples are the methods of Cebeci-Smith 
(ref. 10), Bushnell-Beckwith (ref. 14), Harris (ref. 15), and Herring and 
Mellor (ref. 16). Another approach is to use expressions that consider the rate 
of change of the Reynolds stresses in the governing equations. The methods that 
use this approach are called transport-equation methods. Typical examples are 
the methods of Bradshaw (ref. 17) and Donaldson (ref. 18). For low-speed flows, 
both approaches work equally well. For high-speed flows, however, the mean-flow 
methods seem to be slightly better than the transport-equation methods, chiefly 
because of the inadequate closure assumption accounting for the mean compres- 
sion or dilatation effect. However, a recent report by Bradshaw (Ref. 19) 
seems to substantially improve the predictions of his method for high-speed 
flows. In either case, the governing equations for three-dimensional, com- 
pressible flows are already quite difficult, and there is no need to increase 
the complexity of the equations by using higher-order turbulence models. For 
that reason, we shall restrict our discussion, in this chapter, to the turbu- 
lence models that are based on the eddy-viscosity and mixing-length concepts. 

In particular, we shall describe an eddy-viscosity formulation developed by 
Cebeci (ref. 20), and compare it with others. We shall also present several 
results obtained by that formulation. But first, we shall present a brief 
description of the eddy-viscosity formulation used by Cebeci and Smith for 
two-dimensional flows. 

4.1 Eddy-Viscosity Formulation for Two-Dimensional Compressible Flows 

With Boussinesq's eddy-viscosity concept, we can write the Reynolds shear 
stress, -pu 1 v', as 


-pu 1 V 


T 



(4.1.1) 


*For an excellent discussion of various prediction methods, see a recent 
article by Bradshaw (ref. 13). 
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According to the eddy-viscosity formulation used by Cebeci and Smith in the 
so-called inner region of the boundary layer, e is defined by a modified 
mixing-length expression. In the outer region e is defined by an expres- 
sion based on a velocity defect. For a compressible flow, e is given by 
the following formulas: 


3 U 

ay Y tr 


e - 


a 



Y tr 


o i y i y c 
y c i y i <s 


(4.1.2a) 

(4.1.2b) 


where y c is obtained from the continuity of eddy viscosity. In the above 
equations, L is a modified mixing-length expression given by 

L = <y[l - exp (-y/A)] (4.1.3) 


where 


1/9 1/2 


(4.1.4a) 
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du 

+ 
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e 

dx 
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T 

u 
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1/2 


(4.1.4b) 


(4.1.4c) 


For flows with no mass transfer N can be written as 


N = 


1 - 11.8 — (— ) P + 
y e ' p w' 


|l/2 


(4.1 ,4d) 


According to the study of Cebeci and Mosinskis (ref. 21), the Van Driest 
damping parameter A + and von Karman's parameter k vary with Reynolds number. 
Their variation can be approximated by the following empirical formulas: 

k = 0.40 + ^12 (4.1.5) 


34 


1 + 0.492^ 



14 


(4.1.6) 


A + = 26 + 

_3 

where z 2 e R q x 10 >_ 0.3. 

The parameter a in the outer eddy-viscosity formula is generally 
assumed to be a universal constant equal to 0.0168. According to a recent 
study by Cebeci (ref. 22), however, for values of R. < 6000, a is not a 

universal constant; it varies with R 0 in accordance with the following 

0 

empirical formula: 


1 +n o 

“ ' “o 1 + H 


(4.1.7) 


where a = 0.0168, it = 0.55 and n, which varies from 0 to 1.55 within 
o o 

a R q range of 425 to 6000, is approximated by 


n = 0.55 [1 - exp (-0.243)y 1/2 - 0.298y)], 


y 425 ^ (4.1.8) 


In the definition of y, the R n is defined by 

0 

R .!!& 

\ v w 


where 


e 


k 



(4.1.9a) 


(4.1.9b) 


The low Reynolds number corrections to the eddy-viscosity formulas, 
given by (4.1.5) to (4.1.9), become quite important at high-speed flows. In a 
recent study Bushnell and Morris, (ref. 23), analyzed measurements in hyper- 
sonic turbulent boundary layers at low Reynolds numbers and observed variations 
of the parameters k and a with Reynolds number similar to those given by 
(4.1.5) and (4.1.7). 

The parameter y^ r in the inner and outer eddy-viscosity formulas account: 
for the transitional region that exists between a laminar and turbulent 
boundary layer. It has been used by several authors (refs. 15, 24, 25). 
According to the expression used by Cebeci (ref. 24), it is given by 
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where 


Y tr = 1 - exp 


- 2.68 

3 R a 

U 8 . 

0 _ e tr 

= T ~U~ 

\ A 


- Gr( v>f/ d A(f i 


v tr 


'tr 


(4.1.10) 


A = 60 + 4.68M 


1.92 


(4.1.11a) 


Here, and Re tr are values taken at the start of transition. 


4.2 Extension of the Eddy Viscosity Formulation to Three-Dimensional 

Compres'sib'le Flows 

The eddy-viscosity formulation (4.1.2), which is empirical like all models 
for Reynolds stresses, has worked well for two-dimensional flows. In a recent 
study by Cebeci, Kaups and Mosinskis (ref. 26), it has also been extended to 
handle incompressible three-dimensional flows. Here it will be extended to 
handle compressible flows. In maki-ng this extension, we shall rely heavily 
on our experience with two-dimensional flows and carry over the empirical model 
used for the viscous layer (ref. 21), to three-dimensional compressible flows. 
Because of the lack of data on three-dimensional transitional flows, it is 
difficult to extend the intermittency factor in (4.1.2) to account for the 
transition region. For that reason, the intermittency factor will not be 
included in the formulation of eddy viscosity for three-dimensional flows. 
Furthermore, we shall assume e* = e . 

A Lm 

For the inner region, we shall assume that the inner-eddy-viscosity 
formula is given by the following expression: 


= L 


ft r 

/au\ + /iwY 

\ay/ \3y/ 


il/2 


(4.2.1) 


Here L is given by (4.1.3) and (4.1.4), except that now the friction velocity 
u^ is given by 



where 


w 


w 


= V. 


W 


/au\ 2 + /awY 

\ ay/ lay/. 
x J/ w v J/ ' 


1/2 


(4.2.2) 
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and the dimensionless pressure-gradient parameter p + (uses 2 . 1 . 6 a) and is 
given by 


, vU 3U 

_+ _ s s s 
P “ 3 

u as 

T 


(4.2.3) 


For the outer region, we shall base the eddy-viscosity expression on a 
resultant velocity defect defined by 


u - (u 2 ♦ w 2 )" 2 


and we shall write the outer eddy-viscosity expression as 


e 


0 


= a 


f [u $ - (u 2 + w 2 ) 1/2 ] dy 
0 


(4.2.4) 


Although those inner and outer eddy-viscosity formulas are somewhat 
speculative, they worked quite well for incompressible flow, (ref. 26), and 
are recommended for compressible flows until "better" formulas become available. 

It should be noted that the proposed expressions for the inner and outer 
eddy-viscosity formulas do not, in principle, differ from those suggested by 
Hunt, Bushnell and Beckwith (ref. 27). In a recent study Adams (ref. 28) used 
those transport coefficients in calculating compressible turbulent boundary 
layers on sharp cones at incidence and obtained good agreement with experiment. 


4.3 Attachment-Line Turbulent Flow on an Infinite Swept Wing 

The accuracy of the eddy viscosity presented in Section 4.2 has been 
thoroughly investigated for incompressible infinite swept wings. The calcu- 
lated results agreed well with experiment and with those computed by Bradshaw's 
method (ref. 17). Here, we shall investigate the accuracy of our eddy- 
viscosity formulation for an incompressible attachment-line turbulent flow on 
an infinite swept wing. 

Figure 12 shows a sketch of potential flow streamlines in attachment- line 
region of an infinite swept wing, together with the rectangular coordinate 
system that will be used in this paper. The parameter that determines whether 
the flow will be laminar or turbulent is a dimensionless parameter defined by 
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Figure 12. Sketch of Potential-Flow Streamline in Attachment-Line Region 
of an Infinite Swept Wing and the Coordinate System. 


C* = w. 


dx 


,-1 


(4.3.1) 


It may be regarded as a Reynolds number with the length scale represented by 
the ratio of spanwise velocity, w 


e’ 


to chordwise velocity gradient, du g /dx. 


According to the experiments of Cumpsty and Head (ref. 29), flow along the 
leading edge is fully turbulent for C* >1.4 x 10^. For in5 


the flow is laminar, 
transitional. 


C* < 0.8 x 10' 

In the range 0.8 x 10^ < C* < 1.4 x 10^, the flow is 


4.3.1 Governing Boundary-Layer Equations 

The governing boundary-layer equations for an incompressible turbulent 
flow past a yawed infinite wing, with the use of eddy-viscosity concepts can 
be written as 
Continuity 

&*! r = ° (4 - 3 - 2) 
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Chordwise Momentum 


u 


liL + v — = 

9x ay 


Spanwise Momentum 



(1 + 


+ \ 3U_ 

: ' 9y 


(4.3.3) 


u 



v 


9W _ 9 

3y “ v 3y 


(1 + 


+ n 9W 
: 1 9y 


(4.3.4) 


On the attachment line, u = 0. Therefore, (4.3.3) is singular along the line 
(leading edge) x = 0. To remove the singularity, we differentiate (4.3.3) 

with respect to x and set u and 3v/9x equal to zero. That procedure 
enables (4.3.3) to be written as 


2 9u x 

u + v — - 
x 3y 


1 0.4. 3 

5" t v 

p dx 2 


3y 


+ 3U y 

(1 + * + > 3F 


(4.3.5) 


where u x = 3u/3x. From Bernoulli's equation it follows that at x 


0, 


p dx 7 \ 


du e\ 2 

srl 


(4.3.6) 


Next we introduce a new dependent variable f* defined by 

^du_\-l 




(4.3.7) 


where the prime on f denotes differentiation with respect to the similarity 
parameter n defined by 



(4.3.8) 


with B = (du e /dx) x=Q . 

With (4.3.7) and (4.3.8) we can integrate the continuity equation (4.3.2) 
and can write it as 


v = -/BV f (4.3.9) 

Substituting the expression for v given by (4.3.9) into (4.3.4) and (4.3.5), 
and after performing the necessary transformations, we can write the two 
momentum equations as 
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Chordwise Momentum 


(bf " ) ' + ff" + 1 - (f') 2 = 0 


(4.3.10) 


Spanwise Momentum 

(bg" ) ’ + fg" = 0 (4.3.11) 

In those equations b = 1 + e + and g' denotes the ratio of w/w g . 

Equations (4.3.10) and (4.3.11) are subject to the following boundary 
conditions: 

at n = 0 v 

f = 0 or — - /C*- (mass transfer) 

We - (4.3.12a) 


f' = g = g' = 0 


at n = n 


f = g* = i 


(4.3.12b) 


4.3.2 Eddy-Viscosity Formulation 

The eddy-viscosity formulas (4.2.1) and (4.2.4) become 


e i = (i<y) 2 [1 - exp(-y/A)] 2 


3W 


3y 


£ 0 = “ 


00 

f (w e - w)dy 


0 


(4.3.13a) 


(4.3.13b) 


For zero-pressure gradient flow with no mass transfer, the damping length A 
is 


A = A + v (x w /p )" 1/2 

tn terms of transformed variables (4.3.13) can be written as 


4 = k 2 (C*) 1/2 n 2 ! g" 


1 — exp 


7 w 


1/2 (c .)l/4\l2 


(4.3.14a) 




(4.3.14b) 
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In the study reported here, we have used the above eddy-viscosity formu- 

5 

lation to compute the fully turbulent boundary layers (C* > 1.4 x 10 ) on 
the leading edge of an infinite swept wing. The governing equations, namely, 
(4.3.10) and (4.3.11), were solved by Keller's Box Method. 

When several runs were made for different values of C*, the solutions 
indicated very strong oscillations. The oscillations were small at small 
values of C*, but they became quite strong at high values of C*. It should 
be pointed out that such oscillations are not unusual in turbulent boundary- 
layer calculations. The appearance of such oscillations arise as a result of 
the eddy-viscosity formula given by (4.3.13a); they are observed in all numer- 
ical methods that use (4.3.13a). However, the oscillations in nonsimilar 
turbulent flows are quite small and have no bearing on the accuracy of the 
solutions. 

In order to eliminate the oscillations and provide convergence, we have 
replaced the inner-eddy-viscosity formula (4.3.13a) by another expression, 

e i = <y + [1 - exp (-y/A)]v (4.3.15) 

which, in terms of transformed variables, can be written as 

(4.3.16) 

With that change, no oscillations were observed, and the solutions converged 
quadratically for all values of C* considered. 

4.3.3 Comparison with Experiment 

Detailed measurements of attachment-line flows in turbulent boundary layers 
in incompressible flows are lacking in the literature. The only detailed data 
known to the authors are the data of Cumpsty and Head (ref. 29). For this 
reason, our comparison calculations are limited to that data. Figure 13 shows 
computed and experimental velocity profiles for four values of C*. The 
agreement with experiment is quite satisfactory. 


e i = 
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Figure 13. Comparison of Computed and Experimental Velocity Profiles for the 
Fully Turbulent, Attachment-Line Flow. 


As was mentioned before, the flow is fully turbulent only when 
C* > 1.4 x 10^. For the range of 0.8 x 10^ < C* < 1.4 x 10^, the flow is 
transitional. The calculation for that region was extended by using the inter- 
mittency factor y^ used by Cebeci (ref. 24). For an incompressible flow 
with zero pressure gradient, it is given by 


where 6 is 


v ' 1 - exp 




G = 0.835 x 10' 


-1.34 



(4.3.17) 


(4.3.18) 
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To have similarity we have written (4.3.17) as 


Y tr . = 1 - exp (-Gx* r /w e ) 
which, with the use of (4.3.18), can also be written as 

Y tr = 1 - exp [-0.835 x 10' 3 (R x ) 0 ' 66 ] (4.3.19) 


According to a recent study by Bushnell and Alston (ref. 30), in calculating 

transitional boundary layers it is also necessary to account for the low Reynolds 

number effect (if there is one) in addition to the intermittent behavior of the 

flow. An examination of the experimental data of Cumpsty and Head shows that 

5 5 

for the range of 0.8 x 10 < C* < 1.4 x 10 , The Reynolds number based on e 
varies between 200 and 400. Now the correction to a in (4.1.7), which is for 
a low Reynolds number, applies for R 0 greater than 425. For lower R 0 values, 
we simply extrapolate that curve as shown in figure 14 with a dashed line. The 
resulting (a — R)-curve can be approximated by the following formula: 


a x 10 3 = 194.8 - 128.6 Oog 1() R 0 ) + 30.925 (log in Rj 2 - 2.475 (log in R fi ) 3 


<10V 


for 10 2 < R q < 10 4 . 

0 


1 0 0 ■ 
(4.3.20) 
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10 1 — 
I0 2 


J 1 

I0 3 10* 

R e 


Figure 14. Variation of a with Reynolds Number. 
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Figure 15 shows the transitional boundary layer profiles, together with the 
experimental data of Cumpsty and Head (ref. 29). Those calculations were made by 
multiplying the right-hand side of (4.3.14b) and (4.3.16) by (4.3.19) and by 
varying a in (4.3.14b) as described by (4.3.20). The agreement with experi- 
ment is satisfactory. 

Figure 16 shows a comparison between calculated and measured local skin- 
friction values. Again the agreement with. experiment is satisfactory. 

Finally, we present the computed R Q and H-values in Table 2 at differ- 
ent C*-values. We also present the experimental Revalues given by Cumpsty 
and Head. The agreement between predicted and measured values is quite good. 



Figure 15. Comparison of Computed and Experimental Velocity Profiles for the 
Transitional Attachment-Line Flow. 
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Figure 16. Comparison of Computed and Experimental Skin-Friction Values for 
the Attachment-Line Flow. 


Table 2. 


!. FL and 

H-Values 

for Various 

C*-Vf 

0 

ExR- 

Computed 

x 10“ 5 

R o 

R e 

H 

0.8 

200 

225 

1.76 

1.0 

250 

270 

1.71 

1.2 

295 

313 

1.68 

1.8 

430 

434 

1.60 

2.4 

540 

538 

1.57 

3.0 

640 

634 

1.55 

3.7 

760 

735 

1.53 
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V. OUTLINE OF A GENERAL METHOD FOR COMPUTING COMPRESSIBLE 
THREE-DIMENSIONAL MULTICOMPONENT GAS BOUNDARY LAYERS 


In this chapter, we shall outline a general method for computing compres- 
sible three-dimensional multicomponent gas boundary layers on general configu- 
rations. On the basis of the studies conducted in the earlier chapters, we 
shall give estimates of computation time and computer-storage requirements for 
a typical space-shuttle configuration. Our estimates are given for an equil- 
ibrium or frozen flow consisting of seven-species equations, two momentum 
equations, one energy, and one continuity equation. 

Needless to say, the system of equations under consideration consists of 
highly coupled nonlinear partial-differential equations. They can be solved 
efficiently using the Box Method by following the procedure discussed below. 

In the, discussion,' we shall assume that the governing equations are expressed 
in transformed coordinates. 

1. Express the system in terms of first-order equations. In our case 
this procedure yields 22 first-order equations. 

2. Approximate the system of 22 first-order equations by the difference 
equations for the net in figure 7. 

3. Linearize the resulting nonlinear algebraic equations by Newton's 
method and write them in compound-block-matrix-vector notation as 

A 6. = r • (5.1) 


where 
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The coefficient matrix A is of order 22J + 22 and the vectors 6. 

and r- have this dimension. The blocks A., B., and C. in the 

J J j 

coefficient matrix are of order 22. 

4. Solve the system (5.1) by using the block tridiagonal factorization 
procedure discussed in reference 31. 

To estimate the computer storage and computation time for the Box-scheme 
applied to three-dimensional boundary-layer problems, we suppose there are J 
intervals in the n-di recti on, N intervals in the x-di recti on and K intervals 
in the z-direction. For example, we feel that N = 100, K = 50 and J = 50 
would more than suffice to compute a complete flow field using transformed 
coordinates. The number of basic variables that enter at each net point, 

(x n , z k , rij) Is M h (8 + 2S) where S is the number of species to be 
Included. Specifically, we introduce three variables for each of the x- and 
z-momentum equations, two variables for the energy equation and two variables 
for each species conservation equation (so that each of the equations can be 
reduced to a first-order system). Using S = 7 species yields the M = 22 
basic variable alluded to In steps 1-4 above. 

Since each basic variable requires 4 bytes, it is clear that all of the 
basic variables cannot be stored in the high-speed memory at the same time. 

fi 

This would require for the maximum net cited above, NxKxJx4xM = 22x 10 

3 

bytes or 22 x 10 K-bytes of memory. However, by efficient organization of 

the computer program we need only retain at one time all those basic variables 

on at most 5 "n-columns" [i.e., all points (x , z, , n-) with fixed (x , z.) 

n k j n k 

and all j in 0 <_ j _< J], This requires at most 5 x J x 4 x M = 22 x 10 3 

bytes = 22 K-bytes of high-speed memory. While the solution on one n-column is 

being computed, the data on the next required n-column is being read in from 

auxiliary storage (i.e., disks) and the last computed column is being stored. 

This overlay technique may substantially reduce the delay time in data transfer. 

There is no difficulty in allowing as many as 50 K-bytes for the five columns 

to include fluid properties and other parameters. 

To estimate computation time, we recall that one n-column is obtained by 

solving the linear system (5.1) once for each Newton iteration. The matrix 

elements of A and the inhomogeneous term r. must also be recomputed for 

each iterate. The number of operations to compute these quantities is propor- 
o 

tional to M, J while the number of operations to solve (5.1) is proportional 
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to M J. Thus, we note a very strong dependence on M, the number of basic 
variables to each net point, and only linear dependence on J, the number of 
Intervals across the boundary layer. We recall that the number of n-columns 
is N x K so that the total computation time will also be linear In these 
quantities. Finally, we point out that the computation time is also proportional 
to the average number of Newton iterates employed for each "column." 

For the computations reported in Section 3.2, where M=6, N=25, K = 1 6 
and J = 21 , the total CPU time on an IBM 370/165 was 1.285 minutes. Using 
the above observed linearity with number of n-lntervals if we take J = 50 
rather than J = 21 , the time would be 1.285 x 50/21 = 3 minutes, approximately. 
This Is probably an overestimate since a refinement in the n-spacing would most 
likely reduce the required number of iterations. However, this estimate cor- 

O 

responds to about 7.2 sec/( n ,z)-plane (or 0.9 x 10 sec/net point). Thus, for 
the case of M = 22 (with S = 7 species), N = 25, and K = 16, we estimate: 

7.2 x = 6 mi nutes/(n,z) -plane 

For as many as J = 50 n-points through the boundary layer, this yields five 
hours for the total computation. If we wish to use N = 100 and K = 50 we 
get the tremendous estimate of 125 hours of CPU time. 

This Is unrealistic for practical computations. Fortunately, there are 
several ways In which we can significantly reduce the required computation 
time. First and most basic for Keller's Box-method, is the fact that Richardson's 
extrapolation can be employed and yields two orders of accuracy improvement 
per applications. Thus, with at most N = 50, K = 25 and J = 25 points, 
we can obtain the same accuracy and reduce the computation time by slightly less 
than a factor of 8. So we consider now that 16 hours of CPU time are required 
(with M = 22). 

Another powerful reduction in computation time is obtained by effectively 
reducing the number of basic variables M that are simultaneously coupled in 
solving the Box-difference equations. If only the two momentum and energy 
equations are simultaneously solved, our computing estimate with M = 8 applies. 
Using N = 100, K = 50, J = 50 (i.e., without Richardson's extrapolation), the 

O 

125-hour estimate is now reduced to 125 x (8/22) = 6 hours. However, the 

S = 7 species equations must each be solved and then their updated values 
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employed to recompute the momentum and energy quantities. This "inner-outer" 
iteration procedure should be required at most three times and probably no more 
than twice. Thus, at most 18 hours should be required by this technique. 

If we use the inner-outer iterations only to solve the two momentum equa- 
tions, i.e., M = 6, and then solve separately the energy and species equations, 
the above estimate reduces to about 7.5 hours. This seems to be the most 
promising resolution of the difficulty as now the application of Richardson's 
extrapolation brings us to about one hour of CPU time for an accurate solution. 

We cannot tell at the present stage of development how optimistic or 
pessimistic the above estimates may be. We do feel that they are in the right 
ball park. It clearly should be one of the major objectives of future work in 
this area to test alternatives and to devise the most efficient set of 
procedures. 
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